Many types of physical processes, including fluid flow in a petroleum reservoir, are governed by partial differential equations. These partial differential equations, which can be very complex, are often solved using finite difference, finite volume, or finite element methods. All of these methods divide the physical model into units called gridblocks, cells, or elements. In each of these physical units the solution is given by one or more solution variables or unknowns. Associated with each physical unit is a set of equations governing the behavior of these unknowns, with the number of equations being equal to the number of unknowns. These equations also contain unknowns from neighboring physical units.
Thus, there is a structure to the equations, with the equations for a given physical unit containing unknowns from that physical unit and from its neighbors. This is most conveniently depicted using a combination of nodes and connections, where a node is depicted by a small circle and a connection is depicted by a line between two nodes. The equations at a node contain the unknowns at that node and at the neighboring nodes to which it is connected.
The equations at all nodes are assembled into a single matrix equation. Often the critical task in obtaining the desired solution to the partial differential equations is solving this matrix equation. One of the most effective ways to do this is through the use of incomplete LU factorization or ILU, in which the original matrix is approximately decomposed to the product of two matrices L and U. The matrices L and U are lower triangular and upper triangular and have similar nonzero structures as the lower and upper parts of the original matrix, respectively. With this decomposition, the solution is obtained iteratively by forward and backward substitutions.
The block incomplete LU factorization (BILU) solver has been proved an effective method to solve linear systems resulting from discretization of partial differential equations. In most of the linear systems, the coefficient matrix has dense block sub-matrices as entries. However in the application for oil/gas production surface facilities, those block sub-matrices are sparse. The block sub-matrix might be treated as dense matrix by adding dummy zero terms to the block, but when the block is large and sparse, it might lead to too many dummy entries. For example in a compositional model with 20 hydrocarbon components, such a treatment with dummy entries may increase the matrix size by 10 times, dependent on the hydraulics equations, surface network configurations, active constraint settings and some other factors. Besides that, the number of equations at nodes in some regions may be different from those in other regions, and the resulting block size will vary from region to region. Thus, the conventional BILU with rigid block structures is not a reasonable choice for surface networks.
As shown in FIGS. 1A and 1B, a simple production surface network usually has a tree structure (FIG. 1A) and in some cases includes one or more loops representing, for example, reinjection (FIG. 1B). The governing equations are mass conservation equations on nodes, hydraulics equations on connections, perforation equations on perforation connections and constraint settings on nodes or connections. The primary variables are pressure P and component compositions xi on nodes with two sets of auxiliary variables, total mass flow rate qT on active connections and component flow rates qfi on perforation connections. But the auxiliary variables are eliminated when a Jacobian matrix is assembled, and the final coefficient matrix comprises only mass balance equations with the primary variables P and xi, wherein i=1, . . . , nc−1, and nc is the total number of components. The primary variables are labeled continuously on each node so each node has an unknown vector v consisting of pressure and compositions as follows:
  v  =            [                                                  x              1                                                                          x              2                                                            ⋮                                                              x                              nc                -                1                                                                          P                              ]        .  
The resulting coefficient matrix is illustrated in FIG. 2, which possesses a block structure wherein connection patterns between the blocks represent the physical connections between nodes in the production surface network of FIG. 1B with some virtual connections. The coefficient matrix is not symmetric in incidence. Each row of entries [x] in FIG. 2 represents a mass balance equation and each block row represents a cluster of three mass balance equations related to a respective physical node in FIG. 1B. Because there are seven physical nodes in FIG. 1B, there are seven block rows in FIG. 2. Each entry [x] has a coefficient value and also represents a connection between the physical nodes through the mass balance equation. The elimination of total flow rates may introduce fill-ins in the coefficient matrix as shown in FIG. 2. Because the fill-ins do not exist in the original production surface network connections, they are referred to as virtual connections. For example, in FIG. 2 there is a one way virtual connection from node 5 to node 1, which corresponds to the 15th row in the coefficient matrix, after the total flow rates on connections 1-3 and 3-5 are eliminated although there is no physical connection between node 5 and node 1 in FIG. 1.
Given the coefficient matrix illustrated in FIG. 2, a sparse direct solver may be applied, but it comes with a high cost in computation and storage due to tremendous fill-ins, especially for a matrix resulting from large surface networks with complex configurations and constraint settings. An obvious alternative is the ILU type iterative solver, and its performance may be improved using blocking methods. However, for the production surface network equations, it is difficult to divide the matrix to obtain the block structure with dense or close to dense blocks.